{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "E(Y)=261.53\n"
     ]
    }
   ],
   "source": [
    "#例4-15\n",
    "from scipy.stats import poisson, binom\n",
    "p=poisson.pmf(k=0, mu=1/3)\n",
    "n=365\n",
    "mean=binom.expect(args=(n, p))\n",
    "print('E(Y)=%.2f'%mean)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "E(5X^2)=45.00\n"
     ]
    }
   ],
   "source": [
    "#例4-16\n",
    "from scipy.stats import norm\n",
    "f=lambda x: 5*x**2\n",
    "mean=norm.expect(func=f, loc=0, scale=3)\n",
    "print('E(5X^2)=%.2f'%mean)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "E(XY)=0.125\n"
     ]
    }
   ],
   "source": [
    "#例4-17\n",
    "from scipy.stats import expon\n",
    "meanx=expon.expect(scale=1/2)\n",
    "meany=expon.expect(scale=1/4)\n",
    "print('E(XY)=%.3f'%(meanx*meany))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "E(XY)=0.00\n"
     ]
    }
   ],
   "source": [
    "#练习4-9\n",
    "from scipy.stats import uniform, norm\n",
    "meanx=norm.expect(scale=2)\n",
    "meany=uniform.expect(loc=0, scale=4)\n",
    "print('E(XY)=%.2f'%(meanx*meany))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "E(X)=1.5625\n",
      "D(X)=0.5586\n"
     ]
    }
   ],
   "source": [
    "#例4-18，例4-32\n",
    "import numpy as np\n",
    "from scipy.stats import rv_discrete\n",
    "X=np.array([1,2,3,4])\n",
    "P=np.array([37/64, 19/64,7/64, 1/64])\n",
    "mydist=rv_discrete(values=(X, P))\n",
    "Ex=mydist.expect()\n",
    "Dx=mydist.var()\n",
    "print('E(X)=%.4f'%Ex)\n",
    "print('D(X)=%.4f'%Dx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "E(X)=-0.2\n",
      "E(X^2)=2.8\n",
      "E(3X^2+5)=13.4\n"
     ]
    }
   ],
   "source": [
    "#练习4-10\n",
    "from scipy.stats import rv_discrete\n",
    "import numpy as np\n",
    "X=np.array([-2, 0, 2])\n",
    "P=np.array([0.4,0.3,0.3])\n",
    "mydist=rv_discrete(values=(X, P))\n",
    "dist=mydist()\n",
    "Ex=dist.expect()\n",
    "Ey=dist.expect(func=lambda x:x**2)\n",
    "Ez=dist.expect(func=lambda x:3*x**2+5)\n",
    "print('E(X)=%.1f'%Ex)\n",
    "print('E(X^2)=%.1f'%Ey)\n",
    "print('E(3X^2+5)=%.1f'%Ez)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "E(X)=3/8\n",
      "E(Y)=2\n",
      "E(X^3Y^2)=7/4\n"
     ]
    }
   ],
   "source": [
    "#例4-19\n",
    "import numpy as np\n",
    "from utility import expect, R\n",
    "X=np.array([0, 1])\n",
    "Y=np.array([1, 2, 3])\n",
    "Pxy=np.array([[R(1,4), R(1,8), R(1,4)],\n",
    "              [R(1,8), R(1,8), R(1,8)]])\n",
    "meanx=expect(Pxy, X)\n",
    "g=lambda x, y: y\n",
    "meany=expect(Pxy, Yv=Y, func=g)\n",
    "g=lambda x, y: (x**3)*(y**2)\n",
    "mean=expect(Pxy, X, Y, g)\n",
    "print('E(X)=%s'%meanx)\n",
    "print('E(Y)=%s'%meany)\n",
    "print('E(X^3Y^2)=%s'%mean)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "E(X)=2.00\n",
      "E(Y)=0.00\n",
      "E(Y/X)=-0.07\n"
     ]
    }
   ],
   "source": [
    "#练习4-11\n",
    "import numpy as np\n",
    "from sympy import Rational as R\n",
    "from utility import expect\n",
    "X=np.array([1, 2, 3])\n",
    "Y=np.array([-1, 0, 1])\n",
    "Pxy=np.array([[0.2, 0.1, 0.1],\n",
    "              [0.1, 0.0, 0.1],\n",
    "              [0.0, 0.3, 0.1]])\n",
    "meanx=expect(Pxy, X)\n",
    "g=lambda x, y: y\n",
    "meany=expect(Pxy, Yv=Y, func=g)\n",
    "g=lambda x, y: y/x\n",
    "mean=expect(Pxy, X, Y, g)\n",
    "print('E(X)=%.2f'%meanx)\n",
    "print('E(Y)=%.2f'%meany)\n",
    "print('E(Y/X)=%.2f'%mean)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "E(X)=1500.0,D(X)=375000.0\n"
     ]
    }
   ],
   "source": [
    "#例4-20, 练习4-20\n",
    "from utility import basecont\n",
    "import numpy as np\n",
    "class mydist(basecont):\n",
    "    def _pdf(self, x):\n",
    "        if type(x)!=type(np.array([])):\n",
    "            x=np.array([x])\n",
    "        y=np.zeros(x.size)\n",
    "        d=np.where((x>=0)&(x<1500))\n",
    "        y[d]=x[d]/(1500**2)\n",
    "        d=np.where((x>=1500)&(x<3000))\n",
    "        y[d]=(3000-x[d])/(1500**2)\n",
    "        return y\n",
    "dist=mydist()\n",
    "Ex=dist.expect()\n",
    "Dx=dist.var()\n",
    "print('E(X)=%.1f,D(X)=%.1f'%(Ex,Dx))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "E(X)=0.0\n",
      "D(X)=2.0\n"
     ]
    }
   ],
   "source": [
    "#练习4-12, 例4-33\n",
    "from utility import basecont\n",
    "import numpy as np\n",
    "class mydist(basecont):\n",
    "    def _pdf(self, x):\n",
    "        y=np.exp(-np.abs(x))/2\n",
    "        return y\n",
    "dist=mydist()\n",
    "Ex=dist.expect()\n",
    "print('E(X)=%.1f'%Ex)\n",
    "dist=mydist()\n",
    "Dx=dist.var()\n",
    "print('D(X)=%.1f'%Dx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "E(Y)=0.75\n",
      "E(1/XY)=0.60\n"
     ]
    }
   ],
   "source": [
    "#例4-21\n",
    "from utility import expectcont2\n",
    "f=lambda y, x: 3/(2*x**3*y**2)\\\n",
    "    if (x>1) and (y<x) and (y>1/x)\\\n",
    "    else 0\n",
    "mean=expectcont2(pdf=f, func=lambda y, x: y)\n",
    "print('E(Y)=%.2f'%mean)\n",
    "mean=expectcont2(pdf=f, func=lambda y, x: 1/x/y)\n",
    "print('E(1/XY)=%.2f'%mean)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "E(Y)=1.00\n",
      "E(X)=1.00\n",
      "E(XY)=2.00\n"
     ]
    }
   ],
   "source": [
    "#练习4-13\n",
    "from utility import expectcont2\n",
    "import numpy as np\n",
    "f=lambda y, x: np.exp(-(y+x/y))/y if (x>0)and(y>0) else 0\n",
    "mean=expectcont2(pdf=f, func=lambda y, x: y)\n",
    "print('E(Y)=%.2f'%mean)\n",
    "mean=expectcont2(pdf=f, func=lambda y, x: x)\n",
    "print('E(X)=%.2f'%mean)\n",
    "mean=expectcont2(pdf=f, func=lambda y, x: x*y)\n",
    "print('E(XY)=%.2f'%mean)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "D(X+Y)=5.3333\n",
      "D(2X-3Y)=28.0000\n"
     ]
    }
   ],
   "source": [
    "#例4-31\n",
    "from scipy.stats import norm, uniform\n",
    "Dx=norm.var(scale=2)\n",
    "Ex=norm.expect(scale=2)\n",
    "Dy=uniform.var(scale=4)\n",
    "Ey=uniform.expect(scale=4)\n",
    "print('D(X+Y)=%.4f'%(Dx+Dy))\n",
    "print('D(2X-3Y)=%.4f'%(4*Dx+9*Dy))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "E(X1-X2-2X3+2)=-8.0\n",
      "D(X1-X2-2X3+2)=18.6\n"
     ]
    }
   ],
   "source": [
    "#练习4-17\n",
    "from scipy.stats import norm, binom, poisson\n",
    "Ex1=norm.expect(scale=1)\n",
    "Ex2=binom.expect(args=(10,0.2))\n",
    "Ex3=poisson.expect(args=(4,))\n",
    "print('E(X1-X2-2X3+2)=%.1f'%(Ex1-Ex2-2*Ex3+2))\n",
    "Dx1=norm.var(scale=1)\n",
    "Dx2=binom.var(10, 0.2)\n",
    "Dx3=poisson.var(4)\n",
    "print('D(X1-X2-2X3+2)=%.1f'%(Dx1+Dx2+4*Dx3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "E(A)=8.67\n",
      "D(A)=21.42\n"
     ]
    }
   ],
   "source": [
    "#练习4-18\n",
    "from scipy.stats import uniform\n",
    "Ex=uniform.expect(scale=2)\n",
    "Ex2=uniform.expect(func=lambda x: x**2, scale=2)\n",
    "Ex3=uniform.expect(func=lambda x: x**3, scale=2)\n",
    "Ex4=uniform.expect(func=lambda x: x**4, scale=2)\n",
    "EA=10*Ex-Ex2\n",
    "print('E(A)=%.2f'%EA)\n",
    "EA2=100*Ex2-20*Ex3+Ex4\n",
    "DA=EA2-EA**2\n",
    "print('D(A)=%.2f'%DA)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "E(X)=8.6500, E(Y)=8.6500\n",
      "D(X)=2.5275, D(Y)=2.3275\n"
     ]
    }
   ],
   "source": [
    "#练习4-19\n",
    "from scipy.stats import rv_discrete\n",
    "import numpy as np\n",
    "X=Y=np.array([10,9,8,7,6,5])\n",
    "Px=np.array([0.5,0.1,0.1,0.2,0.05,0.05])\n",
    "Py=np.array([0.4,0.25,0.15,0.05,0.1,0.05])\n",
    "distx=rv_discrete(values=(X,Px))\n",
    "disty=rv_discrete(values=(Y,Py))\n",
    "Ex=distx.expect()\n",
    "Ey=disty.expect()\n",
    "Dx=distx.var()\n",
    "Dy=disty.var()\n",
    "print('E(X)=%.4f, E(Y)=%.4f'%(Ex,Ey))\n",
    "print('D(X)=%.4f, D(Y)=%.4f'%(Dx,Dy))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "D(X)=3/4\n",
      "D(Y)=3/4\n"
     ]
    }
   ],
   "source": [
    "#例4-34\n",
    "import numpy as np\n",
    "from utility import expect\n",
    "from sympy import Rational as R\n",
    "X=np.array([-1, 0, 1])\n",
    "Y=np.array([-1, 0, 1])\n",
    "Pxy=np.array([[R(1,8), R(1,8), R(1,8)],\n",
    "              [R(1,8), R(0), R(1,8)],\n",
    "              [R(1,8), R(1,8), R(1,8)]])\n",
    "Ex=expect(Pxy, X)\n",
    "Ey=expect(Pxy, Yv=Y, func=lambda x, y: y)\n",
    "Ex2=expect(Pxy,X, func=lambda x, y: x*x)\n",
    "Ey2=expect(Pxy, Yv=Y, func=lambda x, y: y*y)\n",
    "Dx=Ex2-Ex**2\n",
    "Dy=Ey2-Ey**2\n",
    "print('D(X)=%s'%Dx)\n",
    "print('D(Y)=%s'%Dy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "D(X)=3/16\n",
      "D(Y)=3/4\n"
     ]
    }
   ],
   "source": [
    "#练习4-21\n",
    "import numpy as np\n",
    "from utility import expect\n",
    "from sympy import Rational as R\n",
    "X=np.array([1, 2])\n",
    "Y=np.array([0, 1, 2, 3])\n",
    "Pxy=np.array([[R(0), R(3,8), R(3,8), R(0)],\n",
    "              [R(1,8), R(0), R(0), R(1,8)]])\n",
    "Ex=expect(Pxy, X)\n",
    "Ey=expect(Pxy, Yv=Y, func=lambda x, y: y)\n",
    "Ex2=expect(Pxy,X, func=lambda x, y: x*x)\n",
    "Ey2=expect(Pxy, Yv=Y, func=lambda x, y: y*y)\n",
    "Dx=Ex2-Ex**2\n",
    "Dy=Ey2-Ey**2\n",
    "print('D(X)=%s'%Dx)\n",
    "print('D(Y)=%s'%Dy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "E(Y)=0.3751\n",
      "D(Y)=0.0593\n"
     ]
    }
   ],
   "source": [
    "#例4-35\n",
    "from utility import expectcont2\n",
    "pdf=lambda y, x: 3*x if(0<x)&(x<1)&(0<y)&(y<x)\\\n",
    "                     else 0\n",
    "Ey=expectcont2(pdf, lambda y, x: y)\n",
    "Ey2=expectcont2(pdf, lambda y, x: y**2)\n",
    "Dy=Ey2-Ey**2\n",
    "print('E(Y)=%.4f'%Ey)\n",
    "print('D(Y)=%.4f'%Dy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "E(X)=0.5556\n",
      "D(X)=0.0802\n"
     ]
    }
   ],
   "source": [
    "#练习4-22\n",
    "from utility import expectcont2\n",
    "pdf=lambda y, x: (x+y)/3 if(0<x)&(x<1)&(0<y)&(y<2)\\\n",
    "                     else 0\n",
    "Ex=expectcont2(pdf, lambda y, x: x)\n",
    "Ex2=expectcont2(pdf, lambda y, x: x**2)\n",
    "Dx=Ex2-Ex**2\n",
    "print('E(X)=%.4f'%Ex)\n",
    "print('D(X)=%.4f'%Dx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "X与Y不相关是True, X与Y相互独立是False\n"
     ]
    }
   ],
   "source": [
    "#例4-37\n",
    "import numpy as np\n",
    "from utility import expect, independent, cov\n",
    "X=np.array([1, 4])\n",
    "Y=np.array([-2, -1, 1, 2])\n",
    "Pxy=np.array([[0, 1/4, 1/4, 0],\n",
    "            [1/4, 0, 0, 1/4]])\n",
    "Ex=expect(Pxy, X)\n",
    "Ey=expect(Pxy, Yv=Y, func=lambda x, y:y)\n",
    "Exy=expect(Pxy, X, Y, func=lambda x, y:x*y)\n",
    "coviar=cov(Exy, Ex, Ey)\n",
    "indep=independent(Pxy)\n",
    "print('X与Y不相关是%s, X与Y相互独立是%s'%(coviar==0, indep))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "X与Y不相关是True, X与Y相互独立是False\n"
     ]
    }
   ],
   "source": [
    "#练习4-24\n",
    "import numpy as np\n",
    "from sympy import Rational as R\n",
    "from utility import expect, independent, cov\n",
    "X=np.array([-1, 0, 1])\n",
    "Y=np.array([-1, 0, 1])\n",
    "Pxy=np.array([[R(1, 8), R(1, 8), R(1, 8)],\n",
    "             [R(1, 8), R(0), R(1, 8)],\n",
    "             [R(1, 8), R(1, 8), R(1, 8)]])\n",
    "Ex=expect(Pxy, X)\n",
    "Ey=expect(Pxy, Yv=Y, func=lambda x, y:y)\n",
    "Exy=expect(Pxy, X, Y, func=lambda x, y:x*y)\n",
    "coviar=cov(Exy, Ex, Ey)\n",
    "indep=independent(Pxy)\n",
    "print('X与Y不相关是%s, X与Y相互独立是%s'%(coviar==0, indep))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Y=0.056*X+2.850\n",
      "X=0.063*Y+2.479\n"
     ]
    }
   ],
   "source": [
    "#例4-40，练习4-26\n",
    "import numpy as np\n",
    "from utility import margDist, expect, rhoxy\n",
    "X=np.array([1, 2, 3, 4, 5])\n",
    "Y=np.array([1, 2, 3, 4, 5])\n",
    "Pxy=np.array([[1/12, 1/24, 1/12, 1/12, 1/24],\n",
    "             [1/24, 1/24, 1/24, 0, 1/24],\n",
    "             [0, 1/24, 1/24, 1/24, 1/24],\n",
    "             [1/24, 1/24, 0, 1/24, 1/24],\n",
    "             [1/30, 1/30, 1/30, 1/30, 1/30]])\n",
    "Ex=expect(Pxy, X)\n",
    "Ex2=expect(Pxy, X, func=lambda x, y: x*x)\n",
    "sigmax=np.sqrt(Ex2-Ex**2)\n",
    "Ey=expect(Pxy,Yv=Y, func=lambda x,y:y)\n",
    "Ey2=expect(Pxy, Yv=Y, func=lambda x, y: y*y)\n",
    "sigmay=np.sqrt(Ey2-Ey**2)\n",
    "Exy=expect(Pxy, X, Y, lambda x, y:x*y)\n",
    "rho=rhoxy(Exy, Ex, Ey, sigmax, sigmay)\n",
    "a=rho*sigmay/sigmax\n",
    "b=Ey-a*Ex\n",
    "print('Y=%.3f*X+%.3f'%(a, b))\n",
    "a=rho*sigmax/sigmay\n",
    "b=Ex-a*Ey\n",
    "print('X=%.3f*Y+%.3f'%(a, b))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD4CAYAAADlwTGnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de5Qc5Xnn8e/T3XO/6DKjC7oYCUlcJO7ICoRs7IM4RhjW8olhLbxkccIx3jX2mpAcFpJdJ8ua2Ng5K7wxOMaAD8GOQZadWHawMQbjNawRSEiC6D5IAkkIjWCkkeau6X72j64e9bR6RiNpqrt6+vc5p89UV1f1vDXT3U/96q1629wdEREpP7FiN0BERIpDBUBEpEypAIiIlCkVABGRMqUCICJSphLFbsDJaG5u9lmzZhW7GSIiJWPt2rXvufukfI+VVAGYNWsWa9asKXYzRERKhpm9NdRjOgQkIlKmVABERMqUCoCISJlSARARKVMqACIiZUoFQESkTKkAiIiUKRUAEZEypQJQJFdccQWPPPJIsZshIlkeeOABPvKRjxS7GQWjAlAk69evZ/PmzcVuhohk2bhxIxs2bCh2MwpGBaBI3J3+/v5iN0NEsvT391NO35KoAlAkqVSKZDJZ7GaISJZkMkkqlSp2MwpGBaBI3F0FQCRiksmkEoCETwlAJHqUAKQglABEokcJQApCncAi0dPf368EIOHK7GEoAYhEixKAhC6zh6ECIBIt6gOQ0CkBiESTEoCETglAJJqUACR0SgAi0aQEIKHL7GHoLCCRaNFZQBI6JQCRaFICkNCpD0AkmtQHIKFTAhCJJiUACZ0KgEg0Zd6T5VIEVACKQJ3AItGUeU+qAEholABEokkJQEKnTmCRaMq8J8ulI1gFoAiUAESiSQlAQqcEIBJNSgASOiUAkWhSApDQ6SwgkWjKvCeVALKY2RIz22pmLWZ2d57Hq8zsqeDx1WY2K+uxe4L5W83smqz5f2ZmG83s38zsB2ZWPRobVAqUAESiSQkgh5nFgQeBa4H5wE1mNj9nsVuBg+4+F1gO3B+sOx9YBiwAlgAPmVnczKYD/xVY6O7nA/FgubKgPgCRaFIfwPEWAS3uvsPd+4AngaU5yywFHg+mVwKLzcyC+U+6e6+77wRagucDSAA1ZpYAaoF3Tm9TSocSgEg0KQEcbzqwO+v+nmBe3mXcvR9oB5qGWtfd9wJ/B7wN7APa3f2X+X65md1mZmvMbM2BAwdG0NzoUwIQiSYlgONZnnm55XGoZfLON7MJpNPBbGAaUGdmN+f75e7+sLsvdPeFkyZNGkFzoy+zd6FOYJFo0VAQx9sDzMy6P4PjD9cMLBMc0hkHtA2z7tXATnc/4O5HgR8Dv38qG1CKlABEokkJ4HivAvPMbLaZVZLurF2Vs8wq4JZg+gbgeU+X0FXAsuAsodnAPOAV0od+Ljez2qCvYDGw+fQ3pzSoD0AkmsqtDyBxogXcvd/MPg88Q/psncfcfaOZ3QuscfdVwKPAE2bWQnrPf1mw7kYzWwFsAvqB2909Caw2s5XAa8H8dcDDo7950aQEIBJN5ZYArJQq3cKFC33NmjXFbsZp27ZtG+eccw719fUcOXKk2M0REdIf+vF4HIB9+/YxderUIrdodJjZWndfmO8xXQlcBLoSWCR6shN5uSQAFYAiUB+ASPRkvx9L6cjI6VABKAIVAJHoUQGQgsjEy1QqVTYvNJGo0yEgKYjsD/1yeaGJRJ0SgBRE9oe+DgOJRIMSgBRE9t6FzgQSiYbs96ISgIRGCUAkepQApCCy9y5UAESiQX0AUhBKACLRowQgBaEEIBI9SgBSENl7F+oEFomG7PeiEoCERglAJHqUAKQg1AcgEj3qA5CCUAIQiR4lACkIJQCR6FECkIJQAhCJHiUAKQidBSQSPToLSApCCUAkepQApCDUByASPeoDkIJQAhCJHiUAKQgVAJHoUQGQglAnsEj0qBNYCkIJQCR6lACkINQJLBI96gSWglACEIkeJQApCCUAkehRApCC0JfCi0SPvhReCkIJQCR6lACkINQHIBI96gOQglACEIkeJQApCCUAkehRApCCUAIQiR4lACkInQUkEj06C0gKQglAJHqUAKQg1AcgEj3qAxiCmS0xs61m1mJmd+d5vMrMngoeX21ms7IeuyeYv9XMrsmaP97MVprZFjPbbGZXjMYGlQIlAJHoUQLIw8ziwIPAtcB84CYzm5+z2K3AQXefCywH7g/WnQ8sAxYAS4CHgucD+AbwC3c/F7gI2Hz6m1MalABEokcJIL9FQIu773D3PuBJYGnOMkuBx4PplcBiM7Ng/pPu3uvuO4EWYJGZNQJ/CDwK4O597n7o9DenNOj7AESiR98HkN90YHfW/T3BvLzLuHs/0A40DbPuWcAB4Ltmts7MHjGzuny/3MxuM7M1ZrbmwIEDI2hu9CkBiESPEkB+lmde7l9nqGWGmp8ALgW+5e6XAJ3AcX0LAO7+sLsvdPeFkyZNGkFzo08FQCR6VADy2wPMzLo/A3hnqGXMLAGMA9qGWXcPsMfdVwfzV5IuCGVBncAi0aNO4PxeBeaZ2WwzqyTdqbsqZ5lVwC3B9A3A854uoauAZcFZQrOBecAr7v4usNvMzgnWWQxsOs1tKRlKACLRU44JIHGiBdy938w+DzwDxIHH3H2jmd0LrHH3VaQ7c58wsxbSe/7LgnU3mtkK0h/u/cDt7p75K38B+H5QVHYAfzLK2xZZSgAi0VOOCeCEBQDA3Z8Gns6Z96Ws6R7gxiHWvQ+4L8/89cDCk2nsWKGhIESiR0NBSEEoAYhETzkmABWAIlAfgEj0lGMfgApAESgBiESPEoAUhBKASPQoAUhBaCgIkejRUBBSEJm9i4qKCiUAkYhIJpNUVFQASgASoszeRSKRUAEQiYjsAqAEIKFRAhCJHiUAKYjM3oUKgEh0JJNJEon0tbFKABIaJQCR6FECkILI7gPQWUAi0dDf368EIOFTAhCJHiUAKQidBSQSPToLSAois3ehAiASHdmdwEoAEhodAhKJHh0CkoLIPg1UncAi0dDf369DQBI+HQISiR4lACmIVCqFmRGPx1UARCJCF4JJQbg7ZkZNTQ3d3d3Fbo6IAN3d3dTW1gJKABKiVCpFLBajvr6ejo6OYjdHRICOjg4aGxsBJQAJUSYB1NfX09nZWezmiAjQ2dlJQ0MDoAQgIcokgLq6OiUAkYjo6OgYKABKABKa7ASgAiBSfMlkku7ubiUACV92H0BnZ2fZ7G2IRFXmUKz6ACR02QkA0JlAIkWWKQBKABK67D4AQIeBRIos8x7M7JQpAUhochOACoBIcWUXADNTApDwZPcBgAqASLHlFgAlAAlNbgLQtQAixZV5D9bX1xOLxZQAJDzqAxCJlsx7sK6uTglAwqU+AJFoyT4EpAQgoVIfgEi0qA9ACkZ9ACLRoj4AKZhMAVAfgEg0dHR0YGZUV1frNFAJV+YQUGbscRUAkeLq6OgYOPwTi8V0CEjCk0kAGhFUJBoyBQBQAshlZkvMbKuZtZjZ3XkerzKzp4LHV5vZrKzH7gnmbzWza3LWi5vZOjP72eluSCnJJABA3wkgEgGdnZ0DBUAJIIuZxYEHgWuB+cBNZjY/Z7FbgYPuPhdYDtwfrDsfWAYsAJYADwXPl/FFYPPpbkSpySQAQENCi0SAEsDQFgEt7r7D3fuAJ4GlOcssBR4PplcCiy39CbcUeNLde919J9ASPB9mNgO4Dnjk9DejtGQnAB0CEim+jo6OgZMylAAGmw7szrq/J5iXdxl37wfagaYTrPsAcBcw7F/azG4zszVmtubAgQMjaG70KQGIRIsSwNAsz7zcv85Qy+Sdb2bXA63uvvZEv9zdH3b3he6+cNKkSSdubQlQH4BItKgPYGh7gJlZ92cA7wy1jJklgHFA2zDrXgl8zMx2kT6kdJWZfe8U2l+SlABEokUJYGivAvPMbLaZVZLu1F2Vs8wq4JZg+gbgeU//BVcBy4KzhGYD84BX3P0ed5/h7rOC53ve3W8ehe0pCeoDEImW7D6AchoKInGiBdy938w+DzwDxIHH3H2jmd0LrHH3VcCjwBNm1kJ6z39ZsO5GM1sBbAL6gdvdPRnStpSM3ARw5MiRIrdIpHy5+3GdwOWSAE5YAADc/Wng6Zx5X8qa7gFuHGLd+4D7hnnuF4AXRtKOsSI7ATQ1NXHw4EGSySTxePwEa4rIaOvo6KCvr4/m5magvBKArgQuguwEMHnyZFKpFG1tbUVulUh5am1tBdLvRSivBKACUATZCWDKlCkA7N+/v5hNEilbmfde5r2oBCChyk4AKgAixZVbAJQAJFRKACLRoQQgBZXbBwAqACLFknnvZS40VQKQUGUngAkTJlBRUaECIFIk+/fvZ+LEiVRUVABKABKy7ARgZkyePHngTAQRKazW1taBwz+gBCAhc/eBBADpY49KACLFsX///kEFQENBSKhSqdRAAgAVAJFiyi0AGgxOQpWbACZPnqwCIFIk+/fvHzgZA5QAJGT5EkBra2vZvOhEoqKnp4fDhw8rAUjhZHcCQ7oA9PX1cejQoSK2SqT85F4DAEoAErLs00ABpk2bBsCePXuK1SSRsrR3717g2HsQlAAkZLkJYO7cuQC8+eabxWqSSFlqaWkBjr0HQQlAQpabAObMmQMcezGKSGG0tLQQi8WYNWvWwDxdCCahyk0AEyZMoKmpSQVApMBaWlo488wzqaysHJinC8EkVLkJANIRVAVApLBaWloGHf4BJQAJWW4CABUAkWLIVwCUACRUQyWAt99+m97e3iK1SqS8tLW1cfDgQSUAKayhEoC7s3PnziK1SqS85DsDCJQAJGRDJQCAbdu2FaNJImVn+/btwLGz8DKUACRU+RLAggULMDPWr19fpFaJlJd169ZRVVXF2WefPWi+EoCEKl8CaGho4Oyzz2bt2rVFapVIeVm7di0XXnjhwBfBZCgBSKjyJQCAyy67TAVApABSqRSvvfYal1122XGPlVMCSBS7AeUolUrx3JYDzLr7XwfNP7yvhoN79zLzC98jXjfhlJ9/11evO90mioxpb775JocPH85bAMyMZDJZhFYVngpAEbg75EkAlVPTHcF977ZQM+eDhW6WSEnJ3YE6GZ2bfgPAl17q4sstg59n/842Fs5sPK22lQodAiqCIQvAlDmA0btve+EbJVJG+t5tgXiCiuYP5HlUg8FJiFKpFMbxBSBWVUvllLPoefv1IrRKpHz0vP06VdPOxeIVxz9oGg5aQjRUAgCoPvMiet/ZQqqvp8CtEikPya52+vbvoHrWxUMuowQgoUnvXQxRAGZdDMl+evdsLGyjRMpEz9tvAE7NmRflX0CngUqYhksAVTPmQzxBzy5dECYShp631mOVNVSecXb+BfSFMBKmVCo1ZAGIVVRTPeN8ut58tWxehCKF4p6iu+VVqs+8CIvF8y5jKAFIiNIf7PkLAEDtuVfS37aHowd2FaxNIuWgd+9mkh3vU3fuHwy9kBKAhCmVSuW9Ejij9uzfB4vRueXFArZKZOzr2vIilqikZs6ioRdSH4CEabg+AIB47TiqP3AhXZv/L+7l8UIUCZunknRteZGasxYSq6odekErn6EgVACKYLg+gIy6CxbTf2gfPW/pmgCR0dC9fTXJzoPUnX/VCZdVAshiZkvMbKuZtZjZ3XkerzKzp4LHV5vZrKzH7gnmbzWza4J5M83s12a22cw2mtkXR2uDSsGJ+gAA6s65klhNI0fWnfrl7iJyzJF1/0q8cdIJh1kxJYBjzCwOPAhcC8wHbjKz+TmL3QocdPe5wHLg/mDd+cAyYAGwBHgoeL5+4M/d/TzgcuD2PM85Zo0kAViikvoLP0L39tUcPfRugVomMjb1HdhFz1sbaLhoyZBn/2RTAjhmEdDi7jvcvQ94Elias8xS4PFgeiWw2NK9nEuBJ9291913Ai3AInff5+6vAbj7EWAzMP30N6c0jCQBADRcdj3EYhx++YfhN0pkDGv/f09hlTXUX3LtiRdWAhhkOrA76/4ejv+wHljG3fuBdqBpJOsGh4suAVbn++VmdpuZrTGzNQcOHBhBc6PvRGcBZSQammm46Bo63vgV/e37C9AykbGn77236dryIg2XXk+8ZgSjfOosoEHyfVLllsehlhl2XTOrB34E3OHuh/P9cnd/2N0XuvvCSZMmjaC50Xeis4CyNV5+IxZLcPDX3w25VSJjj7tz8LnvYFW1NH7w4yNbSdcBDLIHmJl1fwbwzlDLmFkCGAe0DbeumVWQ/vD/vrv/+FQaX6pG0geQkWhopvHyG+ja+iLdb20IuWUiY0t3y2p6dq1j/JWfIl47boRrKQFkexWYZ2azzaySdKfuqpxlVgG3BNM3AM97uoSuApYFZwnNBuYBrwT9A48Cm939f4/GhpSSkfYBZDQu+iMS46fS9ou/J9XXHV7DRMaQZPcR2n75EBXNZ9Jw6ci/Jc+UAI4Jjul/HniGdGftCnffaGb3mtnHgsUeBZrMrAW4E7g7WHcjsALYBPwCuN3dk8CVwB8DV5nZ+uD20VHetshKHwIa+SUYsYoqmj56B/2H9tP2q4fL5sUpcqrcnbZfPkSyq52m6/4Mi5/Elx+WUQEY0V/F3Z8Gns6Z96Ws6R7gxiHWvQ+4L2fei5zMLvAYk+4EPrl1qmeeT+MV/4HDv3uKyiln0XjZvw+ncSJjwJHXfkbXlt8y/g//E1XBV62OmDqBJUwnmwAyxv+7/0jN3N/j4HPfoVvDRYvk1b1jLQef+w41cxfRePkNp/AM5ZMAVACK4FT3LsxiNF//51Q0zeS9f/kKvfu2jXLLREpbz+5/48BPvkrFpFk0X/8X2CnsaKkTWEJ1qgkA0t8bPPmGvyZWXc/+J/87PXs2jXLrREpT9671tK74a+L1Ten3yHADvg2njPoAVACK4HT3LhLjJjPlU/cTr59A64r/Qfebr45Sy0RKU+fWl2hd+T9JTDiDqZ/6ComGplN+LlMfgITJ3U8xmh6TaGxm6k1fJTFhOq0r7+XQSz/Q0NFSdjyV5OAL3+W9f/kKlVPOYspNf0u8bsLpPWkZJYCTODdKRksqlSJ+sqcB5RGvn8DUm79G2zMP0v7i9+ndu4WmJV8YhRaKRN/RQ+/y/s+/Qe/bb1B/8RImLv4slqgYhWcunwSgAlAEJzMUxInEKqppuu5Oqqady8FfP8Y7j36O71zSx6233kospoAnY08ymeTBBx9k32N3gcVo+ugd1F9w9ej9gjJKAPqEKIL03sXoXQZhZjRceh1n/Ok3qZw6l9tuu41Fixbxm9/8ZtR+h0gUPPPMM1x88cV88YtfpGrm+Uy79aHR/fAHXQcg4UongNF/3ooJZzBl2Zd54oknaG1t5cMf/jDXXHMNzz//fNns0cjY4+4899xzXHXVVSxZsoSuri5WrFjB5Bv+hkTj6A8QqaEgJFTpweDC+dObxbj55pvZunUr999/Pxs2bGDx4sUsWrSIFStWcPTo0VB+r8ho6+7u5vHHH+fyyy/n6quvZsuWLSxfvpxNmzZx4403jmhI9VOjBCAhOtnB4E5FTU0Nd911F7t27eLb3/427e3tfPKTn2T69OnceeedvPHGG6H+fpFT4e5s2LCBO+64g+nTp/PpT3+a9vZ2vvWtb7Fjxw7uuOMOqqqqwm1EGSUAdQIXwUi/EOZUzbo793uEp+Mf/zsm7XyNztefZfk3/p7ly5dTMXk2tfOuoPbsK6iYNGvEbdr11ZGPrChyIu7O66+/zsqVK1mxYgXbtm2joqKCT3ziE3z2s5/lQx/6UKjvl+OVTwJQASiwY3sWhR0Lz2Jxaud8kNo5HyTZ1U7nphfo2vIS7S/9gPaX/onE+KlUz7qE6jMvovrMC0f2zUkiwzh+R+SYZFc7PbvW071rHT0715HseB8sRvUHLmDiNZ+n9uwr+F3tOH73i074xdNDPk8olAAkLAMvrILu0QwWrx1H48KlNC5cSrLzIF3bV9PdsprOTS/Qsf7ngFE55SyqP3AhldPOoWraOcQbmgu8FyZjhbuTPNxK794t9O7dTO/ezfTt3wE4sep6qs+8mOrZl1A79/eI140vdnPL6iwgFYACi0IByBavm0DDxUtouHgJnuynd992et5aT89bGzj82k/h1X8eWK5y2jlUTp3LT37Sz/nnn8/s2bN1rYEMkkwm2bFjBxs2bODQb39EX+tO+t7dTrKjDQCrqKLyjLMZ9wefomb2pVROnYvF4kVu9WBWRl8KrwJQYAN7FhEpANksnqB6xnlUzzgPrrwJ7z9KX+sOevdto2/fNnr3baN7+8t8/LffA6C2tpYFCxZwwQUXcN555zFnzpyBW11dXZG3RsJ05MgRWlpaBm7bt29n06ZNvPHGG3R1daUXshgVE2dQfeZFVE07l6rp56b7miL2gZ+PEoCEIrNnYSXwfTiWqKAqOASUkert4uj7u+k78BZHD+zijda3WPvUj0l1Hhq0bqxuPBXjzyAxfiqJcVOJNzQRb2giUZ/+GatpzHtISR3MxefutLW1sWfPHnbv3j3o544dO9i+fTv79+8ftM4ZZ5zBOeecw2c+8xkuuugiLrzwQv7oB29jicoibcVpUB+AhCVqh4BOVqyq9riiAJDq6eDooXfpP/gO/Yfe5ejBffQf2kfP22+QPPICkPOGiieI1x8rCPHaRmI1jXzzmztpamqiubmZ5ubmgena2lMc2ldwd7q6umhvb6etrY0DBw7Q2to65G3v3r10dw/+7ul4PM60adOYPXs2119/PXPnzmXu3LnMmzePOXPmUF9ff9zvtR++W6hNHF0qABKWKB8COh2x6nqqps7N+/V7nuwn2dlG8sj79HekfyY73h+439e6g1TXYVI9R/jCSz/I+/w1NTVMnDiRxsZGGhsbaWhoGJjON6+uro7q6mpqamqorq4eNJ09L5GIzlvA3Ukmkxw9epS+vj56enro7Oykq6tr4JbvfmdnJ4cPH+bQoUMcOnSI9vb246b7+/vz/1KLEatpJF47jnjdOGK106g4/wKqG5qJNzaTaGgm3tBMvG48FovzFvAW8KtDwBpgzR5gT8H+RoWhTmAJSakngFNh8QSJxskkGicz3CU8nkqS6ukg1XWYZHc7qe4jJLsPk+o+TKrrMId6jnCwr5vU/i5Su3fjvd2k+rpJ9XXhvV0clzJGIhbHEpVYLJGejsXAgp+xOFgsfcw6Fk//zNxPtzj44VzygfG4+7C37A/3vr6+gensn6e652kV1cSq6tK36vRPq5pBbOY51M4dPD9eN55Y7XjiteOIVdeXxDH5glInsITl2J5F+RSAkbJYPL0nWjuOCmae1Lrujh/tIdXXjfd2kjrai/cfxZN9eP8wt+RR/GgvpJJ4Kpn+6angfvonnjr2WCoJnsy0mMz/cf27vcGsYJ4F/TyZIT8sPW3xCqw6DnUJLJbA4ulbdSxBTTwoQvEKLBbHKqqwRFX6Z0UVsUH3q4llPaYP8dFjqBNYQjLQCVxGCaAQzAyrrCFWWQP1E4vdHCllZZQAdBJ3gSkBiERc8NYshyKgAlBgx/oA9KcXiaTgvakCIKPu2FlAxW2HiAyvHPoBVAAKTAlAJNpMCUDCUg57FSIlLThBoxzeqyoABXbsLCD96UUiKSgASgAy6sbqlcAiY4cSgISkHPYqREqaEoCE5VgC0J9eJJqUACQk5bBXIVLKTAlAwqIEIBJxOgtIwnLsLKAiN0REhqAEICFRAhCJOCUACcuxvQpFAJFIUh/AYGa2xMy2mlmLmd2d5/EqM3sqeHy1mc3KeuyeYP5WM7tmpM85VpXjF8KIlBJ1AmcxszjwIHAtMB+4yczm5yx2K3DQ3ecCy4H7g3XnA8uABcAS4CEzi4/wOcckXQgmEnXlcwhoJF8IswhocfcdAGb2JLAU2JS1zFLgb4LplcA3LV1GlwJPunsvsNPMWoLnYwTPOWqmTJlCV1dXGE990jIvKg0FIRJRwbernXXWWcRi0XifTp48mTfffHPUn3ckBWA6sDvr/h7g94Zaxt37zawdaArmv5yz7vRg+kTPCYCZ3QbcFtztMLOtI2hzlDQD7+XOfO+nX+e9n369CM0piLzbPMZpm8eYIXYai7LNHR0dp/MtgmcO9cBICkC+35p7cGyoZYaan6+s5j3g5u4PAw8P18AoM7M17r6w2O0oJG1zedA2l76R5Js9MOgbumcA7wy1jJklgHFA2zDrjuQ5RUQkRCMpAK8C88xstplVku7UXZWzzCrglmD6BuB5T3ehrwKWBWcJzQbmAa+M8DlFRCREJzwEFBzT/zzwDBAHHnP3jWZ2L7DG3VcBjwJPBJ28baQ/0AmWW0G6c7cfuN3dkwD5nnP0Ny8SSvbw1WnQNpcHbXOJs3I411VERI4XjXOcRESk4FQARETKlApACMzsL8zMzaw5uG9m9n+CYS9eN7NLs5a9xcy2B7dbhn7WaDKzr5vZlmC7/tnMxmc9NuaHARlL25LNzGaa2a/NbLOZbTSzLwbzJ5rZs8Hr9VkzmxDMH/I1XmqC0QrWmdnPgvuzgyFutgdD3lQG84ccAqdkuLtuo3gjfXrrM8BbQHMw76PAz0lfF3E5sDqYPxHYEfycEExPKPY2nOT2fgRIBNP3A/cH0/OBDUAVMBt4k3SHfzyYPguoDJaZX+ztOMVtHzPbkmfbzgAuDaYbgG3B//RrwN3B/Luz/t95X+OleAPuBP4J+FlwfwWwLJj+B+C/BNOfA/4hmF4GPFXstp/sTQlg9C0H7mLwhW1LgX/0tJeB8WZ2BnAN8Ky7t7n7QeBZ0mMmlQx3/6W79wd3XyZ9TQdkDQPi7juBzDAgA0OLuHsfkBkGpBSNpW0ZxN33uftrwfQRYDPpq/iXAo8Hiz0OfDyYHuo1XlLMbAZwHfBIcN+Aq0gPcQPHb3Pmb7ESWGyncbluMagAjCIz+xiw19035DyUbziN6cPML1V/SnovEMpjm8fStgwpOLRxCbAamOLu+yBdJIDJwWJj5W/xAOkduMxIcE3AoaydnOztGjQEDpAZAqdkjGQoCMliZr8CpuZ56K+AvyR9SOS41fLMG26ojEgZbpvd/SfBMjCCnFkAAAHASURBVH9F+lqP72dWy7P8SQ0DUgJK4v93OsysHvgRcIe7Hx5mB7fk/xZmdj3Q6u5rzezDmdl5Fh3uSz1KaptVAE6Su1+db76ZXUD6WPeG4E0yA3jNzBYx/JAYH86Z/8KoN/o0DbXNGUHn9fXAYg8OiDL8cB9jZRiQMT2kiZlVkP7w/767/ziYvd/MznD3fcEhntZg/lj4W1wJfMzMPgpUA42kE8F4M0sEe/nZ25XZ5j05Q+CUjmJ3QozVG7CLY53A1zG4g+yVYP5EYCfpDuAJwfTEYrf9JLdzCekrvSflzF/A4E7gHaQ7TRPB9GyOdZwuKPZ2nOK2j5ltybNtBvwj8EDO/K8zuBP4a8F03td4qd5I75hlOoF/yOBO4M8F07czuBN4RbHbfbI3JYDCeJr0WRItQBfwJwDu3mZm/4v02EgA97p7ae1BwDdJf8g/GySfl939P3sZDAPiQwyTUuRmjZYrgT8G3jCz9cG8vwS+Cqwws1uBt4Ebg8fyvsbHiP8GPGlmXwbWkR76BoYYAqeUaCgIEZEypbOARETKlAqAiEiZUgEQESlTKgAiImVKBUBEpEypAIiIlCkVABGRMvX/AUItpI1NkTgQAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "E(V)=-0.0000, D(V)=72200.0000\n"
     ]
    }
   ],
   "source": [
    "#程序4.19，4.20\n",
    "from scipy.stats import rv_continuous\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt#导入绘图对象plt\n",
    "class vdist(rv_continuous):\n",
    "    def _cdf(self, v):\n",
    "        y=np.zeros(v.size)\n",
    "        d=np.where((v>-380)&(v<380))\n",
    "        y[d]=np.arcsin(v[d]/380)/np.pi+1/2\n",
    "        d=np.where(v>=380)\n",
    "        y[d]=1\n",
    "        return y\n",
    "dist=vdist()\n",
    "v=np.linspace(-500,500, 256)\n",
    "plt.plot(v, dist.pdf(v), color='black')\n",
    "V=dist.rvs(size=10000)\n",
    "plt.hist(V, density=True, histtype='stepfilled')\n",
    "plt.show()\n",
    "Ev=dist.expect()\n",
    "Dv=dist.var()\n",
    "print('E(V)=%.4f, D(V)=%.4f'%(Ev, Dv))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD4CAYAAADlwTGnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXhUVbb38e/KTBLCEKKgoIBMRkCGgCitNoIMomB3qwzyChqgRVDBRhvae23boS94HVHURrBFWiZpB1qQwQnEAQgOQUAwIGIUIRKEDGRe7x8pvDEmpCDDrmF9nofHqlP7nPpVWalVZ+9z9hFVxRhjTPAJcR3AGGOMG1YAjDEmSFkBMMaYIGUFwBhjgpQVAGOMCVJhrgOcjCZNmmjLli1dxzDGGL+yZcuWH1U1ofxyvyoALVu2JCUlxXUMY4zxKyLyTUXLrQvIGGOClBUAY4wJUlYAjDEmSFkBMMaYIGUFwBhjgpQVAGOMCVJWAIwxJkhZATCmBhUVFXHs2DFsmnXjD6wAGFMNhw8fZu7cuQwZMoQ2bdoQGRlJdHQ0UVFRdOnShcmTJ7N69WpKSkpcRzXmV8SffqkkJSWpnQlsfMGBAwe49957ef755ykoKKB169YkJSXRrl07YmJiOHz4MFu2bOGDDz4gLy+PDh06MHXqVMaMGUNoaKjr+CbIiMgWVU0qv9yvpoIwxrWSkhIef/xx/vrXv5KXl0dycjJjx46le/fuiMiv2ufn5/Pvf/+b//3f/2Xs2LHMmTOHefPm0bFjRwfpjfkl6wIyxkuHDh1iyJAh/OlPf+LSSy9l27ZtPPvssyQlJVX45Q8QGRnJyJEj+eSTT1i4cCF79uyhW7duzJ49u47TG/NrVgCM8cKuXbvo3r07a9eu5amnnuI///kP7dq183p9EWHEiBFs376d/v37M2nSJG699VaKiopqMbUxJ2ZdQMZUITU1lcsvv5ySkhI2bNhAjx49TnlbCQkJvP766/z5z3/mkUceIT09naVLlxIeHl6DiY3xju0BGHMCW7du5dJLLyUiIoL333+/Wl/+x4WGhvLwww/zxBNP8NprrzFq1CjbEzBO2B6AMZXYt28fAwcOJCYmhvfff5+avhjRbbfdRkFBAXfeeSfR0dE8//zzlY4lGFMbrAAYU4HDhw8zaNAgsrOz2bBhQ41/+R83depUsrOz+dvf/kaHDh3485//XCvPY0xFvCoAIjIQeAIIBeaq6oxyj0cCLwLdgUPAMFXdKyLxwDKgB/CCqk6qYNvLgdaqasfFmWprOW1FtbehWkLGsvs4tncXp193H1e9tA/Yd8J19s4YfMrP99e//pWdO3cyffp0EhMTueqqq055W8acjCrHAEQkFJgNDAISgREikliuWTJwWFXbAI8BMz3L84D/BqZWsu3fA9mnFt2Y2nH0o5c5tieFxn3HEXV251p/PhFh3rx5dOvWjZEjR/LVV1/V+nMaA94NAvcE0lR1j6oWAIuBoeXaDAXme24vA/qKiKhqjqpuoLQQ/IKIxAJ3AA+ccnpjatixbz7npw0vEZ14KbFdr6iz542OjubVV18lPDyckSNHUlBQUGfPbYKXNwXgTODbMvfTPcsqbKOqRcARIL6K7d4PPALkepXUmFpWnJfNoTceJazRGcQPmFTnA7ItWrRg7ty5pKSkcM8999Tpc5vg5E0BqOivoPwEQt60+b/GIl2ANqr6apVPLjJeRFJEJCUjI6Oq5sacssNv/YPinMM0uWoqIRH1nGT4/e9/z/jx43nooYd47733nGQwwcObApAOtChzvznwfWVtRCQMaABknmCbFwLdRWQvsAFoJyLvVdRQVeeoapKqJiUkJHgR15iTl7vrI3K2vUuDC4cR2bSN0yyPPfYYrVu3Zty4cRw7dsxpFhPYvCkAm4G2ItJKRCKA4cDycm2WA6M9t68B3tETTDOqqs+o6hmq2hL4DbBLVX97suGNqQkledkcWjObiNPPocFF17mOQ3R0NHPmzCEtLY17773XdRwTwKosAJ4+/UnAamAHsFRVt4nIfSIyxNNsHhAvImmUDuxOO76+51f+o8AYEUmv4AgiY5z66f0FlOQepfHAW5FQ35iS4bLLLiM5OZlHHnmETz75xHUcE6DsegAmoJzseQD5+7/ihxfvoH73K2nc74+n/LzVOQ+gMocPH+bcc8/l7LPP5qOPPiIkxGZuMaemsusB2CfKBC0tKSZzzdOExjSk4cWjXMf5lUaNGvHQQw+xadMmXnzxRddxTACyAmCCVs629yj44SsaXZZMSGSM6zgVGjVqFL169WLatGkcPXrUdRwTYKwAmKBUUpjPT+8vIKJZW6LPvdR1nEqFhIQwa9YsDh48yP333+86jgkwVgBMUMpKeZ3irB9p1CfZ52fg7NGjB2PGjGHWrFns23fiOYmMORlWAEzQKc49wpGPX6ZemwuIauEfcxD+7W9/Q0TssFBTo6wAmKBz5INFaGE+jX47xnUUr7Vo0YKJEycyf/58tm/f7jqOCRBWAExQKcz8jqzP3iT2/AGEx7eoegUfMn36dGJiYviv//ov11FMgLACYILKkQ8WIaFhNPzNSNdRTlqTJk248847efXVV9m0aZPrOCYAWAEwQaMw8ztydqynftfBhMY0ch3nlEyePJmEhASmTZuGP53EaXyTXRLSBI0jHy5GwsKJ6/n7Gt92TVyJzFtFna/m3befo9n1M4g6q3OtnIVsgoPtAZigUJj5HTnb11G/yxWExjR0HadaYs8fSEhMQ458uNR1FOPnrACYoHDko6VIaDhxF9T8r/+6FhIeSVyP35H3zWfkf7/TdRzjx6wAmIBXeHg/OdveJbbLQL/t+y+vfpdBhETFcuSjJa6jGD9mBcAEvNJf/2E0uOAa11FqTEhkNPWThnIsbROpqamu4xg/ZQXABLSirB9Lf/137k9obGD8+j+ufverkIh6/P3vf3cdxfgpKwAmoGVt+Q9oCfV7XO06So0LjYqlftfBLF26lJ07bSzAnDwrACZgleTnkvXpm0S37014w6au49SKuB5XExkZycMPP+w6ivFDVgBMwMr+fBVakFsrx/37itCYhtxwww0sWLCAjIwM13GMn7ECYAKSFhdxNGU5kWd1IrJZW9dxatWUKVPIz8/n6aefdh3F+BmvCoCIDBSRnSKSJiLTKng8UkSWeB7fKCItPcvjReRdEckWkafKtI8WkRUi8qWIbBORGTX1gowByNmxnuKsH2kQwL/+j+vQoQODBw9m9uzZ5OXluY5j/EiVBUBEQoHZwCAgERghIonlmiUDh1W1DfAYMNOzPA/4b2BqBZt+WFU7AF2B3iIy6NRegjG/pKoc3fQK4U3OIqr1r66DHZDuuOMOMjIyeOmll1xHMX7Emz2AnkCaqu5R1QJgMTC0XJuhwHzP7WVAXxERVc1R1Q2UFoKfqWquqr7ruV0AfAI0r8brMOZned98TmHGXuJ6/M7nr/ZVU/r06cP555/Po48+apPEGa95UwDOBL4tcz/ds6zCNqpaBBwB4r0JICINgauAtyt5fLyIpIhIig1yGW9kffIGIfXiiEn03Wv91jQR4Y477mD79u2sWbPGdRzjJ7wpABX9hCr/E8ObNr/esEgYsAiYpap7KmqjqnNUNUlVkxISEqoMa4Jb0ZGDHEvbROz5A5CwCNdx6tTw4cNp1qwZjz76qOsoxk94UwDSgbKXTmoOfF9ZG8+XegMg04ttzwG+UtXHvWhrTJWyPi2dlrl+1+AbUoqIiOCWW25hzZo17Nq1y3Uc4we8KQCbgbYi0kpEIoDhwPJybZYDoz23rwHe0So6IkXkAUoLxeSTi2xMxY4dO0b252uIbtuLsLjTXMdxYuzYsYSHh9shocYrVRYAT5/+JGA1sANYqqrbROQ+ERniaTYPiBeRNOAO4OdDRUVkL/AoMEZE0kUkUUSaA3dTelTRJyLymYiMrckXZoLPokWLKMnLon73K11HcaZp06b84Q9/4IUXXiAnJ8d1HOPjvLoimKquBFaWW3ZPmdt5wLWVrNuyks0Gx+EZQagur451nKqy/4UHCW9yNpEtOtX58/uSiRMnsnjxYhYuXMi4ceNcxzE+zM4ENgEh/7sdFB7cQ/3uVwbNoZ+V6d27N507d2b27Nl2SKg5ISsAJiBkbfkPIZExxCT2cR3FORFh4sSJfP7553z44Yeu4xgfZgXA+L2irEPk7vqQmM6XExIR5TqOT7j++uuJi4tj9uzZrqMYH2YFwPi97M9XQUkJ9bsOdh3FZ8TExDBmzBiWLVvGgQMHXMcxPsoKgPFrWlJMdupbRLXsQnijZq7j+JRbbrmFwsJC5s6d6zqK8VFWAIxfy9v7GcVZGcSeP8B1FJ/Tvn17+vXrx7PPPktRUZHrOMYHWQEwfi07dS0h9eKIbnOB6yg+aeLEiaSnp7NiRd0fmmt8nxUA47eKc4+Q+9XHxJzXBwkLdx3HJ1155ZU0a9aM5557znUU44OsABi/lbPtXSgpIrbz5a6j+KywsDBuvPFG3nzzTdLT013HMT7GCoDxS6pKduoaIpq1IyKhpes4Pi05OZmSkhKef/5511GMj7ECYPxSwf5dFP64j9jO/V1H8XmtW7emX79+zJs3j+LiYtdxjA+xAmD8UnbqGiQ8kphzL3EdxS+MGzeOffv2sXbtWtdRjA+xAmD8TklBHjk71hPd/jeEREa7juMXhg4dSpMmTWww2PyCFQDjd3J3bkALjtng70mIjIxk9OjRLF++3M4MNj+zAmD8TnbqGsIan0lk8/NcR/ErY8eOpaioiBdeeMF1FOMjrAAYv1J4KJ389O3Edro86Kd9PlkdOnTg4osvZu7cuTZNtAGsABg/k711LUgIsR37uo7il8aNG0daWhrvvfee6yjGB1gBMH5Di4vI/uJt6rXpSWhsI9dx/NI111xDw4YNbTDYAFYAjB85tieFkpyfbPC3GurVq8eoUaP497//zaFDh1zHMY55VQBEZKCI7BSRNBGZVsHjkSKyxPP4RhFp6VkeLyLviki2iDxVbp3uIrLVs84ssQ5dU4Xs1DWExjamXusk11H82rhx4ygoKGDBggWuoxjHqiwAIhIKzAYGAYnACBFJLNcsGTisqm2Ax4CZnuV5wH8DUyvY9DPAeKCt59/AU3kBJjgUZR3i2O4UYjpehoSEuo7j1zp37kzPnj157rnnbDA4yHmzB9ATSFPVPapaACwGhpZrMxSY77m9DOgrIqKqOaq6gdJC8DMRaQbEqepHWvoJfBG4ujovxAS2nC/eBi0htpN1/9SEsWPHsn37djZu3Og6inHImwJwJvBtmfvpnmUVtlHVIuAIEF/FNstOTVjRNgEQkfEikiIiKRkZGV7ENYFGVcneupbIFh0Jb1zhx8ScpGHDhhEdHc28efNcRzEOeVMAKuqbL7/f6E2bU2qvqnNUNUlVkxISEk6wSROo8r/9gqLD+23wtwbFxcVx3XXXsXjxYrKzs13HMY54UwDSgRZl7jcHvq+sjYiEAQ2AzCq22byKbRoDeCZ+i4gmun1v11ECSnJyMtnZ2bz88suuoxhHvCkAm4G2ItJKRCKA4cDycm2WA6M9t68B3tETjC6p6n4gS0R6eY7+uQF4/aTTm4BXkp9D7s4PiUm8hJDwKNdxAkrv3r1p3769dQMFsSoLgKdPfxKwGtgBLFXVbSJyn4gM8TSbB8SLSBpwB/DzoaIishd4FBgjIulljiCaAMwF0oDdwJs185JMIMnZvg4tyrd5/2uBiHDTTTfxwQcf8OWXX7qOYxzw6jwAVV2pqu1U9RxVfdCz7B5VXe65naeq16pqG1Xtqap7yqzbUlUbq2qsqjZX1e2e5Smq2tGzzUkn2mMwwSs7dS3hCS2JaNrWdZSAdMMNNxAaGmpXCwtSdiaw8VkFB7+m4IeviO1sE7/VlqZNm3LllVcyf/58CgsLXccxdcwKgPFZ2alrIDSMmPP6uI4S0JKTkzl48CArVqxwHcXUMSsAxidpUSE5294juu2FhNaLcx0noA0aNIhmzZrZYHAQsgJgfFLuVx9Rkpdlx/7XgbCwMEaPHs3KlSv5/ns7GjuYWAEwPik7dS2hcacR1bKL6yhB4aabbqKkpIT58+dX3dgEDCsAxucUHTlA3t5Pie3UDxH7iNaFtm3bcskll/D888/bBHFBxP66jM/JTn0LEGI793MdJagkJyeTlpbG+vXrXUcxdcQKgPEpWlJM9ta3iGrZhbC401zHCSrXXHMNcXFxdk5AELECYHxK3t7PKM7KsDN/HYiOjmbEiBG8/PLLHDlyxHUcUwesABifkp26lpB6cUS37eU6SlBKTk7m2LFjLF682HUUUwesABifUZx7hNyvPiYm8bdIWLjrOEEpKSmJTp062TkBQcIKgPEZOdveg5IiYs+37h9XRITk5GQ2b97M1q1bXccxtcwKgPEJqkp26hoimrUlIqGl6zhBbdSoUURERNheQBCwAmB8QsH+XRT++I0N/vqA+Ph4rr76ahYsWEB+fr7rOKYWWQEwPiE7dS0SFknMuZe4jmIoHQzOzMzk9dftOk2BzAqAca6kII+cHeuI7tCbkMgY13EM0LdvX1q0aGHdQAHOCoBxLnfnB2jBMev+8SGhoaHceOONrF27lm+++cZ1HFNLrAAY57JTVxPW6Awim5/nOoop48YbbwTghRdecBvE1Jow1wFMcCs8lE5++nYaXjrGrvp1ilpOq70LuUSedT73P/o0z+d2+9XEfHtnDK615zV1w6s9ABEZKCI7RSRNRKZV8HikiCzxPL5RRFqWeWy6Z/lOERlQZvkUEdkmIl+IyCIRiaqJF2T8S/bWtSAhxHbs6zqKqUBs58spPnqQvL2fu45iakGVBUBEQoHZwCAgERghIonlmiUDh1W1DfAYMNOzbiIwHDgPGAg8LSKhInImcBuQpKodgVBPOxNEtLiI7C/ept45PQiNbeQ6jqlAdLsLCYmKLb08pwk43uwB9ATSVHWPqhYAi4Gh5doMBY5fSWIZ0FdK9+eHAotVNV9VvwbSPNuD0u6neiISBkQDdimiIHNsTwolOT/Z4K8Pk7AIYs7rQ+5XH1F87KjrOKaGeVMAzgS+LXM/3bOswjaqWgQcAeIrW1dVvwMeBvYB+4EjqlrhTwwRGS8iKSKSkpGR4UVc4y+yU9cQGtOIeuckuY5iTiC28+VQXFQ6VYcJKN4UgIpG5spfMqiyNhUuF5FGlO4dtALOAGJEZFRFT66qc1Q1SVWTEhISvIhr/EFRdibHdqcQ07EvEhLqOo45gYjTWhPRtA3ZqWvsamEBxpsCkA60KHO/Ob/urvm5jadLpwGQeYJ1+wFfq2qGqhYCrwAXncoLMP4p54u3QUvsou9+IrZzfwoz9lLwQ5rrKKYGeVMANgNtRaSViERQOli7vFyb5cBoz+1rgHe09KfCcmC45yihVkBbYBOlXT+9RCTaM1bQF9hR/Zdj/IFqCdmfryGyRUfCG5fvTTS+KObcS5CwCBsMDjBVFgBPn/4kYDWlX9JLVXWbiNwnIkM8zeYB8SKSBtwBTPOsuw1YCmwHVgETVbVYVTdSOlj8CbDVk2NOjb4y47Py9m2l6Kf9xJ4/oOrGxieERMUS3b43OdvXUVKY5zqOqSFenQimqiuBleWW3VPmdh5wbSXrPgg8WMHyvwJ/PZmwJjBkf76GkMgYottZr58/ie18OTnb3iV354fEdrzMdRxTA2wqCFOnio8dJXfXB8Sc14eQ8EjXccxJiGzRibCGzawbKIBYATB1KueLd6G4yLp//JCIENv5cvK//YLCw3baTiCwAmDqjKqS/flqIpq1I+K0Vq7jmFMQ0/EykBCyU9e6jmJqgBUAU2fyv/uSwkP77MxfPxZWvwn1Wncn54u3KSoqch3HVJMVAFNnslNXI+FRdtUvPxfb+XKKszNZtWqV6yimmqwAmDpRkp9L7pfvE3PuJYRERruOY6qh3jk9CYluaFcLCwBWAEydyNn+HlqYb4O/AUBCw4jteBlvvPEGBw4ccB3HVIMVAFMnsj9fTXhCSyKatXMdxdSA2M6XU1RUxIsvvug6iqkGKwCm1uX/kEbBgd3Enj/ArvoVIMLjW3DRRRcxb948myDOj1kBMLUuO3XNz/PKm8CRnJzMzp07+eijj1xHMafICoCpVSUFeeRse4/o9r0JjYp1HcfUoOuuu47Y2FgbDPZjVgBMrcrZsQ4tyLXB3wAUGxvLsGHDWLJkCVlZWa7jmFNgBcDUGlUl+9OVhDc5m8jm57mOY2pBcnIyOTk5LF261HUUcwqsAJhaU7B/FwUHdlO/6xU2+BugevXqRWJiIs8995zrKOYUWAEwtSbr05VIRD0b/A1gIsIf//hHNm7cyCeffOI6jjlJVgBMrSg+drT0zN/E39qZvwHuhhtuIDo6mtmzZ7uOYk6SFQBTK3K2voUWFVC/6xWuo5ha1rBhQ0aNGsXChQvJzMx0HcecBCsApsaplpD12ZtEnplo0z4HiYkTJ5KXl8cLL7zgOoo5CVYATI3L2/s5RYf3U7+b/foPFp07d+Y3v/kNTz/9NCUlJa7jGC95VQBEZKCI7BSRNBGZVsHjkSKyxPP4RhFpWeax6Z7lO0VkQJnlDUVkmYh8KSI7ROTCmnhBxr2sT1cQEt2A6Ha9XUcxdeiWW25h9+7drFljl4z0F1UWABEJBWYDg4BEYISIJJZrlgwcVtU2wGPATM+6icBw4DxgIPC0Z3sATwCrVLUDcD6wo/ovx7iWnp7OsbRNxHa+HAkLdx3H1KE//OEPnH766TYY7Ee82QPoCaSp6h5VLQAWA0PLtRkKzPfcXgb0ldIDv4cCi1U1X1W/BtKAniISB1wCzANQ1QJV/an6L8e49o9//ANUqd9lkOsopo5FREQwbtw4VqxYwddff+06jvGCNwXgTODbMvfTPcsqbKOqRcARIP4E67YGMoB/isinIjJXRGIqenIRGS8iKSKSkpGR4UVc40peXh7/+Mc/qNemB2ENTncdxzgwfvx4RIRnn33WdRTjBW8KQEWncJaf/7WyNpUtDwO6Ac+oalcgB/jV2AKAqs5R1SRVTUpISPAirnFl0aJFZGRkUD+p/A6iCRYtWrRg6NChzJs3j7y8PNdxTBW8KQDpQIsy95sD31fWRkTCgAZA5gnWTQfSVXWjZ/kySguC8VOqyuOPP06nTp2IOquz6zjGoUmTJnHo0CEWLVrkOoqpgjcFYDPQVkRaiUgEpYO6y8u1WQ6M9ty+BnhHS68SsRwY7jlKqBXQFtikqj8A34pIe886fYHt1XwtxqF169aRmprK7bffbvP+BLk+ffrQqVMnHnvsMbtYjI+rsgB4+vQnAaspPVJnqapuE5H7RGSIp9k8IF5E0oA78HTnqOo2YCmlX+6rgImqWuxZ51bgJRFJBboAf6+5l2Xq2uOPP06TJk0YOXKk6yjGMRFhypQpbN26lbffftt1HHMC4k8VOikpSVNSUlzHMOXs2bOHNm3acPfdd3P//ffTctoK15FMHdg7Y3Clj+Xl5XH22WeTlJTEihX2eXBNRLaoalL55XYmsKm2J598ktDQUCZMmOA6ivERUVFR3HLLLaxcuZIvv/zSdRxTCSsAplqOHj3KvHnzGDZsGGeccYbrOMaHTJgwgcjISB5//HHXUUwlrACYannhhRfIysri9ttvdx3F+JjTTjuNUaNG8eKLL3Lo0CHXcUwFrACYU1ZcXMysWbO46KKL6NGjh+s4xgdNnjyZY8eOlZ4hbnyOFQBzyl577TV2797NlClTXEcxPqpjx47079+fp556ioKCAtdxTDlWAMwpUVVmzpzJOeecw+9+9zvXcYwPmzJlCvv372fJkiWuo5hy7DDQAFabh2Pm7dvKgUXTadz/FrvqV5A60WGgZakqHTt2JCwsjM8++8xOFHTADgM1NerIxmWERDckpmNf11GMjxMR7rzzTlJTU3nzzTddxzFlWAEwJ63g4Nfk7dlCXPerCAmPdB3H+IGRI0fSokUL/ud//sd1FFOGFQBz0o5uegUJjyLWun6MlyIiIpg6dSobNmxgw4YNruMYDysA5qQUHT1Izo71xJ4/gNB69V3HMX4kOTmZ+Ph4ZsyY4TqK8bACYE7K0Y3/BiCuh835b05OTEwMt99+OytWrCA1NdV1HIMVAHMSirIOkfX5GmI79iUs7jTXcYwfmjRpErGxscycOdN1FIMVAHMSjm56BUqKibvwOtdRjJ9q1KgRf/zjH1m8eDFfffWV6zhBzwqA8UpxzmGyP1tFzHl9CG/Y1HUc48emTp1KREQEDz74oOsoQc8KgPHK0U2vosWFNLBf/6aamjZtyoQJE/jXv/5FWlqa6zhBzQqAqVJx7hGyPl1J9LkXE974TNdxTAC46667iIiI4P7773cdJahZATBVOrr5NbQwnwYXDnMdxQSIsnsBNhbgjhUAc0LFOYfJ2vIfojv8hogmZ7mOYwLInXfeSWRkJA888IDrKEHLqwIgIgNFZKeIpInItAoejxSRJZ7HN4pIyzKPTfcs3ykiA8qtFyoin4rIG9V9IaZ2HPnoZbSogIYXj3IdxQSYsnsBu3btch0nKFVZAEQkFJgNDAISgREikliuWTJwWFXbAI8BMz3rJgLDgfOAgcDTnu0ddzuwo7ovwtSOoqMHyfpsJbGd+lnfv6kVd911F1FRUdxzzz2uowQlb/YAegJpqrpHVQuAxUD500CHAvM9t5cBfaV0ztehwGJVzVfVr4E0z/YQkebAYGBu9V+GqQ1HPlgMQIPewx0nMYHq9NNPZ8qUKSxZsoQtW7a4jhN0vCkAZwLflrmf7llWYRtVLQKOAPFVrPs4cBdQcqInF5HxIpIiIikZGRlexDU1ofBQOtlb36J+18F21q+pVXfeeSfx8fFMnz7ddZSg400BqOjqDeWvIlNZmwqXi8iVwEFVrbLkq+ocVU1S1aSEhISq05oa8dOGl5CwCBr0utZ1FBPgGjRowN13383atWt56623XMcJKt4UgHSgRZn7zYHvK2sjImFAAyDzBOv2BoaIyF5Ku5QuE5F/nUJ+Uwvy9+8i98v3iUsaSmhMQ9dxTBCYMGECZ511FtOmTaOk5ISdAqYGeVMANgNtRaeEzagAABI/SURBVKSViERQOqi7vFyb5cBoz+1rgHe09FqTy4HhnqOEWgFtgU2qOl1Vm6tqS8/23lFVO8zEB6gqh9+ZS0h0Q+Iu+IPrOCZIREVFcf/997NlyxaWLVvmOk7QqLIAePr0JwGrKT1iZ6mqbhOR+0RkiKfZPCBeRNKAO4BpnnW3AUuB7cAqYKKqFtf8yzA1JXfXh+Snb6fhxaMIiYx2HccEkeuvv55OnToxffp08vLyXMcJCl6dB6CqK1W1naqeo6oPepbdo6rLPbfzVPVaVW2jqj1VdU+ZdR/0rNdeVX91QVBVfU9Vr6ypF2ROnRYV8tN7/yQ8oSWxnS93HccEmdDQUB599FH27NnD448/7jpOULAzgc3Pjm75D0U//UCjPslISGjVKxhTw/r168fQoUN54IEH+P778kONpqZZATBA6YRvRz5aQr3WSdRr1dV1HBPEHn74YQoLC/nLX/7iOkrAswJgAPhp3Xy04BiN+iS7jmKCXJs2bZgyZQrz589n06ZNruMENCsAhvzvdpCduoa4HlcT3qRF1SsYU8vuvvtumjZtym233WaHhdYiKwBBTkuKObTmaULrN6FB7xGu4xgDQP369ZkxYwYbN25k3rx5ruMELCsAQS7rkzcoPPg1jfuOJySinus4xvzshhtu4NJLL+Wuu+7iwIEDruMEJCsAQawo6xA/vf8volp3p167C13HMeYXRIRnn32W3NxcpkyZ4jpOQLICEMQOv/0cWlxE4343Uzp5qzG+pUOHDkyfPp1FixaxevVq13ECjhWAIJXz5QZyd26gYe8RhDdq5jqOMZWaPn067dq1Y8KECeTm5rqOE1DCXAcwda849wiZa58homkbm+/HnLKW01bU2XPldR/DgUV/oelvr+foplfr7HkDne0BBKHMtc9SkpdD/BWT7Yxf4xeizupMbNcryNr8OuvWrXMdJ2BYAQgyOTs/IPfL92nQezgRCS1dxzHGa41+exNhDZsyZswYsrKyXMcJCFYAgkhxzmEy1zxDxOnn0OCCa1zHMeakhEREET94Cvv27eOOO+5wHScgWAEIEqol/LjicUryc4gfPAUJteEf43+imidy5513MnfuXFasqLsxiEBlBSBIZG1+nbyvt9D4srHW9WP82t/+9jc6derEjTfeaDOGVpMVgCCQ/0Mah9fNp167C4nteoXrOMZUS2RkJIsXLyYnJ4eRI0dSVFTkOpLfsgIQ4Eryc/lx+UxCYxoSP/A2O+HLBITExESeeeYZ1q1bx7333us6jt+yAhDAVJVDKx+n6KcDNLlqKqH16ruOZEyNueGGG7jpppv4+9//bmcJnyIrAAHs6Mcvk7vrQxr9dgxRLTq6jmNMjXvyySc577zzGDVqFN98843rOH7HqwIgIgNFZKeIpInItAoejxSRJZ7HN4pIyzKPTfcs3ykiAzzLWojIuyKyQ0S2icjtNfWCTKmVK1fy0/oFRCdeSv0ev3Mdx5haER0dzbJlyygsLGTIkCFkZ2e7juRXqiwAIhIKzAYGAYnACBFJLNcsGTisqm2Ax4CZnnUTgeHAecBA4GnP9oqAP6nquUAvYGIF2zSnaNeuXYwcOZKI01sTP/BW6/c3Aa19+/YsXbqUL774glGjRtkFZE6CN3sAPYE0Vd2jqgXAYmBouTZDgfme28uAvlL6rTMUWKyq+ar6NZAG9FTV/ar6CYCqZgE7gDOr/3LMjz/+yFVXXUV4eDgJv7ubkPAo15GMqXX9+/fnscce4/XXX+fuu+92HcdveFMAzgS+LXM/nV9/Wf/cRlWLgCNAvDfrerqLugIbK3pyERkvIikikpKRkeFF3OCVk5PDlVdeyb59+3jttdcIa3Ca60jG1Jlbb72V8ePHM2PGDP75z3+6juMXvCkAFfUfqJdtTriuiMQC/wYmq+rRip5cVeeoapKqJiUkJHgRNzgVFRUxfPhwNm/ezKJFi+jdu7frSMbUKRHhqaee4vLLL2fs2LG89tprriP5PG8KQDpQ9krhzYHyp9/93EZEwoAGQOaJ1hWRcEq//F9S1VdOJbwppapMmDCBN954g9mzZ3P11Ve7jmSME+Hh4bzyyiv07NmTYcOG8e6777qO5NO8KQCbgbYi0kpEIigd1F1ers1yYLTn9jXAO6qqnuXDPUcJtQLaAps84wPzgB2q+mhNvJBgparcfvvtzJ07l7vvvpubb77ZdSRjnIqNjWXFihW0bduWIUOGkJKS4jqSz6qyAHj69CcBqykdrF2qqttE5D4RGeJpNg+IF5E04A5gmmfdbcBSYDuwCpioqsVAb+D/AZeJyGeefzZHwUlSVSZPnsyTTz7JlClTuP/++11HMsYnNG7cmDVr1tCkSRP69+/P5s2bXUfySVL6Q90/JCUlqVXzUse//GfNmsWUKVN45JFHfnW4Z11escmYurJ3xmDv2+7dy2WXXcahQ4dYuXJl0I6NicgWVU0qv9zOBPZDRUVF3HzzzcyaNYvJkydX+OVvjIGWLVuyfv16mjZtyoABA2xMoBybFL6W1fSv8JLCPH5c/hDH0jYR1+taXonoy6vTV9bocxgTSJo3b866devo168fV1xxBQsWLOCaa+yCSGB7AH6lOPcIBxbdzbG0zTS+/GYaXTrafvkb44WmTZvy3nvv0bVrV6699lpmzpyJP3V/1xYrAH6i4ODX/LBgKgUH95Dwu+nU73al60jG+JUmTZrwzjvvMGzYMKZNm8a4ceMoLCx0Hcsp6wLyAznb3+PQm08SEhVD0xF/J/LMc11HMsYvRUVFsXDhQtq0acODDz7Ijh07WLJkCc2bN3cdzQnbA/BhWlxI5tvP8eN/Hiai6Tk0G/2EffkbU00hISE88MADLFy4kM8//5wuXbqwatUq17GcsALgowoy9rL/xTvISnmd+t2v4vThfyc0tpHrWMYEjBEjRpCSkkKzZs0YNGgQ06dPJz8/33WsOmXnAdSykz0KSLWErM2vcXj9i4RExhI/8Fai215QS+mMMSWFeRx+aw7ZqWsIb3IW8VdMIbJZ2wrbnsw5CL6ksvMAbAzAh+T/kEbmmmco2L+Tem17ET9gEqExDV3HMiaghYRHET/oNuq1u5DMVU/yw4I/EXfB72lw0fCAn07dCoAPKMnL5qf3/0XWpysJqRdH/JV/Iibxt3aIpzF1KPqcHkQlP03mO/M4+vEycrato1GfG4nucHHA/i1aAXCopDCf7E9XcOTjZZTkZVO/6xU0vHgUIVGxrqMZE5RComJpcsXtxHbqS+Zbc/hx+UNEfrKCRn1uIvKM9q7j1TgrAA5oUSHZX7zNkQ8WUZx9iKiWXWn02zFEnH6O62jGGCCqRUeajX6M7NS1/PT+An5Y8CfqtU5i0+8T6Nmzp+t4NcYKQB0qycsm67NVZG1ZTnF2JhFntKfJVX8i6qzOrqMZY8qRkFDqdxlIzLmXkPXpCo5ufIULLriAAQMGMHnyZPr3709IiH8fSGkFoA7k/5BGduoacra9ixYcI+rsLsQPup2oVt0Ctm/RmEAREhlNg17XUr/rYCaclsasWbMYNGgQ7du359Zbb+X666+nYUP/PFjDDgOtJfv372fZsmVMffAJCg7sRsIiiO7wG+KSribi9Nau4xljTsHeGYMpKChg6dKlPPHEE6SkpBAZGcnQoUMZPXo0/fv3JyzM935XV3YYqBWAGrR7927efPNNXn75Zd5//31UlfDTWlP//AHEJF5qg7vG+Lmy5wGoKlu2bGH+/PksXLiQzMxMGjVqxODBgxk6dCgDBgygfv36DtP+HysAteCnn35i/fr1rF69mtWrV7N7924AzjvvPK699lquvfZarnjxa8cpjTE1pbITwQoKCli5ciWvvvoqb7zxBpmZmURERNCnTx/69OnDpZdeSvfu3QkPD6/jxKWCvgBUd15+LS6iMDOd/O93kv/dlxR8v5PCQ/sAkPAoos7qRFSrbtRr1Y3wxmdW67mMMb7JmzOBi4qK+OCDD3j99ddZtWoVO3bsACAmJoaLLrqIXr160a1bN7p160aLFi3qZBzQCsBJFoD873eS983nFGTspTDjGwozv4OSIgBC6sUReUZ7Is5oT9SZiUSeeS4S5qayG2PqzqlMBXHgwAHWr1/PunXrWL9+Pdu2baOkpASA+Ph4unTpQmJiIomJiVx33XU0bty4pmNXrwCIyEDgCSAUmKuqM8o9Hgm8CHQHDgHDVHWv57HpQDJQDNymqqu92WZF6rIAHF73Akc/XkZo3GlEJJxNeEJLIhLOJqJpW8IanWFH7xhjTklJYR6FGd9QcGA3BT+kUXDwawoPfYsW5nHmhH8SFpfwq3WqOwfRKc8FJCKhwGzgciAd2Cwiy1V1e5lmycBhVW0jIsOBmcAwEUkEhgPnAWcAb4lIO886VW3Tqbiev6dBr+sIiYx2HcUYE0BCwqOIPKP9L84sVlWKszIIrd+kbrN40aYnkKaqe1S1AFgMDC3XZigw33N7GdBXSn8iDwUWq2q+qn4NpHm25802nQqtF2df/saYOiEihMWdVuc9C94csHom8G2Z++lA+fmJf26jqkUicgSI9yz/uNy6x0dIq9omACIyHhjvuZstIju9yFyXmgA/ug5xEvwtL/hfZn/LC/6X2d/yQjUyy8xqP/fZFS30pgBUVJLKDxxU1qay5RXteVQ4GKGqc4A5JwrokoikVNS35qv8LS/4X2Z/ywv+l9nf8oJvZvamCygdaFHmfnPg+8raiEgY0ADIPMG63mzTGGNMLfKmAGwG2opIKxGJoHRQd3m5NsuB0Z7b1wDvaOnhRcuB4SISKSKtgLbAJi+3aYwxphZV2QXk6dOfBKym9JDN51V1m4jcB6So6nJgHrBARNIo/eU/3LPuNhFZCmwHioCJqloMUNE2a/7l1Qmf7Z6qhL/lBf/L7G95wf8y+1te8MHMfnUimDHGmJrj35NZG2OMOWVWAIwxJkhZAThFIrJERD7z/NsrIp95lrcUkWNlHnvWdVYAEblXRL4rk+uKMo9NF5E0EdkpIgNc5jxORP5XRL4UkVQReVVEGnqW++T7e5yIDPS8j2kiMs11nvJEpIWIvCsiO0Rkm4jc7lle6efDF3j+xrZ6sqV4ljUWkbUi8pXnv41c5wQQkfZl3sfPROSoiEz2xffYxgBqgIg8AhxR1ftEpCXwhqp2dJvql0TkXiBbVR8utzwRWETp2dlnAG8B7Y4P1rsiIv0pPZqsSKT0NBhV/bOvvr/w87QpuygzxQkwwpemOBGRZkAzVf1EROoDW4Crgeuo4PPhK0RkL5Ckqj+WWfYQkKmqMzzFtpGq/tlVxop4PhPfUXqi64342HtsewDV5Jny4jpKv0T9UWXTdTilqmtUtchz92NKzxXxdT4/xYmq7lfVTzy3s4Ad/N/Z+f6m7BQ08yktZL6mL7BbVb9xHaQiVgCq72LggKp+VWZZKxH5VETWicjFroJVYJKnS+X5MrvLFU314WtfCDcBb5a576vvrz+8lz/z7E11BTZ6FlX0+fAVCqwRkS2e6WEATlfV/VBa2IDTnKWr3HB++ePQp95jKwAnICJvicgXFfwr+6tuBL/8H7wfOEtVuwJ3AAtFJM4H8j4DnAN08WR85PhqFWyqTvoFvXl/ReRuSs8hecmzyNn76wVn7+XJEpFY4N/AZFU9SuWfD1/RW1W7AYOAiSJyietAVZHSk1yHAC97Fvnce+x7Vy/2Iara70SPS+m0F7+n9DoIx9fJB/I9t7eIyG6gHVDr17KsKu9xIvIc8IbnrrNpObx4f0cDVwJ9PWeWO31/veAXU5yISDilX/4vqeorAKp6oMzjZT8fPkFVv/f896CIvEppd9sBEWmmqvs9YxsHnYb8tUHAJ8ffW198j20PoHr6AV+qavrxBSKS4Bn4QURaUzr9xR5H+X7m+QM57nfAF57blU3X4ZSUXjDoz8AQVc0ts9wn318Pn5/ixDNmNQ/YoaqPllle2efDORGJ8QxYIyIxQH9K85WdgmY08LqbhJX6Re+AL77HtgdQPeX79wAuAe4TkSJKr4J2s6pm1nmyX3tIRLpQ2iWxF/gjnHi6DseeAiKBtaXfWXysqjfju+9vpdOmOI5VXm/g/wFbxXPoMvAXYERFnw8fcTrwqudzEAYsVNVVIrIZWCoiycA+4FqHGX9BRKIpPRqs7PtY4d+gS3YYqDHGBCnrAjLGmCBlBcAYY4KUFQBjjAlSVgCMMSZIWQEwxpggZQXAGGOClBUAY4wJUv8fIUpd5TLMQekAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#程序4.21\n",
    "from scipy.stats import norm \n",
    "x=np.array([np.random.choice(a=V, size=100).mean()\n",
    "            for k in range(500)])\n",
    "plt.hist(x, density=True, histtype='stepfilled')\n",
    "v=np.linspace(-80, 80,256)\n",
    "plt.plot(v, norm.pdf(v, scale=np.sqrt(Dv)/10), color='black')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
